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Abstract 

Semiclassical  approximations  for  quantum  time  correlation  functions  are 
presented  for  both  electronically  adiabatic  and  nonadiabatic  dynamics  along 
with  discussions  of  the  operator  ordering  and  the  classical  limit.  With  the 
combined  use  of  the  initial-value  representation  of  the  semiclassical  propaga¬ 
tor,  a  discrete  algorithm  to  evaluate  the  Jacobi  matrices,  semiclassical  op¬ 
erator  ordering  rules,  and  the  stationary-phase  filter  technique,  a  practical 
algorithm  is  developed  to  calculate  quantum  time  correlation  functions.  This 
approach  holds  considerable  promise  for  simulating  the  quantum  dynamics  of 
realistic  many-body  systems.  Some  simple  illustrative  examples  are  used  to 
demonstrate  the  feasibility  and  accuracy  of  the  algorithm. 
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I.  INTRODUCTION 


It  is  well-known  from  linear  response  theory  that  the  response  of  a  system  to  a  weak 
external  force  can  be  formulated  in  terms  of  a  time  correlation  function  for  the  relevant 
dynamical  variable  of  the  system  (see,  e.g.,  Refs.  1-3).  Therefore,  time  correlation  functions 
play  a  central  role  in  the  study  of  dynamical  processes,  such  as  chemical  reactions,  light 
scattering  spectra,  spectroscopic  lineshapes,  transport  properties,  etc.  Classically,  the  evo¬ 
lution  of  the  system  obeys  the  Newtonian  equation  of  motion,  which  serves  as  the  basis  for 
molecular  dynamics  (MD)  simulations.  Quantum  mechanically,  the  probabilistic  wavefunc- 
tion  propagates  according  to  the  Schrodinger  equation,  which  in  principle  cannot  be  solved 
by  means  of  deterministic  trajectory  dynamics.  Due  to  the  importance  of  time  correlation 
functions,  much  effort  has  been  devoted  to  the  development  of  methods  to  calculate  them 
quantum  mechanically;  unfortunately,  few  methods  have  been  successful  in  applications  to 
realistic  many-body  quantum  systems.  In  fact,  it  turns  out  that  real-time  quantum  propaga¬ 
tion  is  a  truly  formidable  numerical  problem  because  large  sign  fluctuations  in  the  real-time 
propagator  overwhelm  the  contribution  from  the  physical  quantities  of  interest  (see,  e.g.. 
Refs.  4-10).  Thus,  to  this  day  the  real-time  propagation  of  many-body  quantum  systems 
remains  a  daunting  challenge. 

There  have  been  several  attempts  to  calculate  quantum  time  correlation  functions  ex¬ 
actly  using  the  Feynman  path  integral  formulation.'*'®’^^  For  example,  by  virtue  of  the  nu¬ 
merical  matrix  multiplication  method  (NMM),*^  Thirumalai  and  Berne  were  able  to  cal¬ 
culate  the  symmetrized  dipole-dipole  time  correlation  function  for  a  proton  moving  in  a 
one-dimensional  bistable  potential.*^  While  the  NMM  approach  becomes  prohibitive  for 
many-dimensional  systems,  it  is  also  fruitless  to  directly  apply  Monte  Carlo  methods  to 
evaluate  time  correlation  functions  in  such  systems  due  to  large  phase  cancellations.  To 
treat  the  generic  problem  of  performing  many-dimensional  averages  of  highly  oscillatory 
integrands— which  are  the  origin  of  the  difficulty  in  direct  Monte  Carlo  calculations  of  such 
functions-several  stationary  phase  Monte  Carlo  (SPMC)  methods  have  been  developed. 
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The  implementation  of  these  and  other  techniques  makes  it  possible  in  some  cases  to  sim¬ 
ulate  the  dynamics  of  simplified  many-dimensional  quantum  systems.  For  example,  there 
have  been  a  series  of  studies  on  electron  transfer  dynamics  as  represented  by  the  spin-boson 
model  and  its  multi-state  generalization. For  system-bath-type  Hamiltonians  having 
harmonic  baths,  quasi-adiabatic  propagator  path  integral  methods  (QUAPI)  have  also  been 
developed  which  propagate  adiabatically  a  one-dimensional  system  in  which  the  harmonic 
bath  has  been  incorporated  through  an  analytic  influence  functional.^^  By  virtue  of  this 
algorithm  and  discrete  variable  representation  (DVR)  quadrature,  a  detailed  study  of  quan¬ 
tum  rates  for  a  double  well  coupled  to  a  harmonic  bath  was  recently  presented^®  along 
with  a  comparison  to  approximate  theories.  Unfortunately,  all  of  the  methods  described 
above  are  either  not  directly  applicable  to  “real”  nonlinear  many-body  potentials  or  become 
numerically  intractable  for  anything  but  the  short  time  dynamics  of  such  systems. 

One  “exact”  alternative  to  the  direct  real  time  quantum  dynamics  approach  is  based  on 
the  fact  that  real-time  correlation  functions  can  be  formally  related  to  their  imaginary-time 
counterparts  through  analytic  continuation  {l/ksT  =  /?  — ^  Thus,  in  principle, 

one  can  simulate  a  quantum  system  with  an  equilibrium  path  integral  Monte  Carlo  method 
at  several  values  of  imaginary  time  and  infer  the  real-time  quantities  through  the  analytic 
continuation  (see,  e.g..  Refs.  31—34).  In  practice,  however,  the  analytic  continuation  is 
rather  sensitive  to  statistical  fluctuations  in  imaginary  time  data  so  this  approach  has  suf¬ 
fered  from  numerical  instabilities.  Gubernatis  and  coworkers  have  recently  introduced  the 
maximum  entropy  method  (MEM)  which  appears  to  improve  the  stability  of  the  analytic 
continuation.^'^®'^®  The  MEM  has  proven  to  be  reliable  and  efficient  in  similar  ill-posed  in¬ 
version  problems,  so  its  application  in  path  integral  simulations  seems  novel  and  promising. 
Using  this  technique,  Gallicchio  and  Berne^^  have,  for  example,  calculated  the  dipole  absorp¬ 
tion  spectrum  of  an  electron  in  fluid  helium  and  found  good  agreement  with  some  previous 
analytic  results.  The  implementation  of  this  and  other  versions  of  the  MEM  allows  one  to 
evaluate  the  lower  frequency  portion  of  the  absorption  spectra  with  good  accuracy,  but  it 
probably  requires  further  effort  to  determine  the  high-frequency  portion  which  is  essential  in 


describing  short-to-intermediate  time  quantum  dynamics,  e.g.,  photodissociation  processes, 
optical  control  experiments,  quantum  tunneling,  and  charge  transfer. 

As  an  alternative  to  the  numerical  evaluation  of  the  exact  quantum  time  propagation  in 
many-body  systems,  one  can  develop  approximate  methods  for  quantum  dynamics  on  which 
stable  and  feasible  numerical  algorithms  can  be  based  to  compute  time  correlation  func¬ 
tions.  One  such  approach^®“^^  has  been  developed  by  the  present  authors  and  is  based  on 
the  dynamical  properties  of  the  centroid  variable  in  Feynman  path  mtegration.^’^^  In  this  ap¬ 
proach,  called  “Centroid  Molecular  Dynamics”  (CMD),  a  quasiclassical  dynamics  algorithm 
is  employed  to  compute  an  approximation  to  the  the  Kubo-transformed  quantum  dynamical 
time  correlation  function.  There  are  now  several  encouraging  results  from  applications  of 
CMD  to  a  variety  of  non-trivial  systems.'‘i'^'  The  simplicity  and  stability  of  this  method 
makes  CMD  a  promising  candidate  for  quantum  dynamical  simulations  in  the  condensed 

phase  where  other  methods  become  impractical. 

In  the  present  paper,  however,  a  different  and  promising  approach  for  the  calculation 
of  approximate  quantum  dynamical  time  correlation  functions  will  be  developed  based  on 
semiclassical  arguments,  some  of  which  originate  from  the  earliest  formulations  of  the  “old 
quantum  theory.  Indeed,  since  those  early  days  many  attempts  have  been  made  to  elucidate 
and  utilize  the  relationship  between  classical  dynamics  and  its  quantum  counterpart.  In 
time-independent  quantum  mechanics,  this  is  commonly  known  as  the  WKB  (Wentzel- 
Kramers-Brillouin)  approximation  for  one-dimensional  problems,  and  it  can  be  generalized 
to  many-dimensional  problems  as  in  classical  S-matrix  theory'*®  '**  (Miller-Marcus  theory) 
and  EBK  (Einstein-Brillouin-Keller)  quantization  theory  (see,  e.g..  Refs.  3,  49).  On  the 
other  hand,  time- dependent  semiclassical  mechanics  was  first  studied  by  Van  Vleck“  and 

later  extended  by  many  others.'*®  '*®’®* 

Although  semiclassical  approaches  have  found  wide  use  in  various  analytical  theories,  the 
development  of  semiclassical  quantum  dynamics  as  a  numerical  algorithm  has  been  hindered 
by  two  major  drawbacks:  the  root  search  problem  and  the  caustics  problem.  These  two  dif¬ 
ficulties  can  be  avoided  in  some  cases,  e.g.,  in  the  context  of  Miller  s  S-matrix  theory. 
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with  the  help  of  an  initial-value  representation  in  which  an  integration  in  phase  space  by¬ 
passes  the  root  search.  Moreover,  Campolieti  and  Brumer  have  extended  the  initial-value 
analysis  to  real-time  propagation  and  have  thus  suggested  a  semiclassical  approach  in  which 
the  classical  trajectories  evolve  according  to  the  initial  phase  space  representation.  Earlier, 
Miller  and  Heller  proposed  an  initial-value  propagation  of  wavefunctions  which  introduces 
integrations  over  initial  and  final  positions  and  thus  allows  for  a  change  of  variables  to  the 
initial  phase  variables. 

In  the  present  paper,  we  show  how  the  initial-value  semiclassical  approach  for  comput¬ 
ing  the  quantum  propagator  can  be  used  to  calculate  time  correlation  functions.  To  be 
more  specific,  we  have  re-derived  the  initial-value  semiclassical  propagator  from  a  discrete 
perspective  and  found  an  alternative  for  evaluating  the  Jacobi  matrices  in  the  discretized 
formalism.  Importantly,  these  new  developments  allow  us  to  formulate  the  theory  and  a 
tractable  numerical  algorithm  for  both  adiabatic  and  nonadiabatic  semiclassical  time  prop¬ 
agation  of  the  nuclei  in  quantum  systems.  This,  along  with  an  initial-value  expression  for 
the  evaluation  of  quantum  operators,  makes  it  possible  to  implement  semiclassical  dynamics 
in  the  calculation  of  quantum  dynamical  time  correlation  functions.  The  emphasis  in  the 
present  paper  is  on  a  formulation  amenable  to  realistic  many-body  simulations. 

The  sections  of  this  paper  are  organized  as  follows:  In  Sec.  II,  the  semiclassical  approx¬ 
imation  for  quantum  time  correlation  functions  is  described  and  rederived  in  the  adiabatic 
dynamics  limit  from  both  the  boundary- value  and  initial-value  perspectives,  the  latter  being 
shown  to  be  superior  for  our  purposes.  This  derivation  is  next  generalized  in  Sec.  Ill  to  the 
nonadiabatic  limit.  Then,  in  Sec.  IV  a  stationary-phase  filter  method  is  introduced  to  aid  in 
the  actual  implementation  of  the  initial-value  semiclassical  methodology  and  some  numerical 
examples  are  studied  in  Sec.  V  to  demonstrate  the  feasibility  of  the  algorithm.  Concluding 
remarks  are  given  in  Sec.  VI,  while  the  Appendices  contain  important  supporting  material. 
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II.  SEMICLASSICAL  THEORY:  ADIABATIC  DYNAMICS 


A.  Boundary-value  Formulation 


1.  Van  Vleck  Formula  for  the  Propagator 


It  is  well-known  that  semiclassical  mechanics  can  be  largely  understood  as  an  asymp¬ 
totic  analysis  of  functional  integrals  in  terms  of  which  to  second  order  is  equivalent 
to  the  stationary  phase  approximation.'^®  In  this  subsection,  the  boundary-value  Van  Vleck 
formulation  of  time-dependent  semiclassical  theory  will  be  reviewed  for  completeness  and  as 
background  material  for  subsequent  developments.  To  start,  the  real-time  propagator  can 
be  expressed  according  to  Feynman’s  prescription  of  path  integrals  as 

G{qi,cit]t)  =  (qt|e'''^*/^|qi)  =  j  Vq{t')  exp{iS[q{t')]/h}  ,  (2.1) 

where  the  action  S[q{t')],  given  by 

S[q{t')]  =  dt'L[q{t'),q{t')]  =  j\dt'  {^q(0  '  “  '  “  '^"[^(0]}  >  (2-2) 

is  evaluated  with  the  Lagrangian  L[q{t'),  q(t')]  along  the  path  q(t)  subject  to  the  boundary 
conditions 


q(0)  =  qi 

q{t)  =  qt  .  (2-3) 


Following  common  notation,  fonts  with  hats  denote  operators  and  bold  fonts  denote  vectors 
or  matrices;  in  particular,  the  vectors  q  =  {^i,  92,  •  ■  •  >  and  p  =  {pi,P2,  -  ■  ■  ,Pn}  represent, 
respectively,  V-dimensional  coordinates  and  their  conjugate  momenta  in  an  iV-degree-of- 
freedom  system,  whereas  m  is  the  diagonal  mass  matrix.  An  application  of  the  stationary- 


phase  approximation  to  Eq.  (2.1)  gives® 


Gsciqu  q<;  0  -  E  j  ^  det  ^ 


ex 


(2.4) 
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where  the  summation  is  carried  over  all  possible  stationary  paths,  and  Sst  =  qti  i)  is 

the  classical  action  associated  with  a  given  stationary  path.  The  stationary  phase  condition 
[6S/6q{t%t  =  0  determines  the  classical  trajectory,  thus  leading  to  the  Euler-Lagrange 


equations 


d  fdL\  dL  ..  ,  dV 


(2.5) 


with  the  boundary  conditions  in  Eq.  (2.3). 

In  the  short  time  limit,  the  determinant  in  Eq.  (2.4)  is  positive  and  this  semiclassical 
expression  is  exactly  the  original  Van  Vleck  short-time  propagator.  In  general,  the  deter¬ 
minant,  termed  the  Van  Vleck  determinant,  can  be  written  in  a  more  useful  form  as 

=  \detJq{t)\~^^'^ex^{-iTTu{t)/2)  ,  (2.6) 

where  u{t),  known  as  the  Maslov  index,^^  is  the  number  of  sign  changes  of  the  determinant  as 
the  trajectory  evolves  in  time  from  0  to  t.  The  Jacobi  matrices,  defined  as  3q{t )  —  5q(t  )/5pi 
and  Jp(t')  =  dp{t')/dpi,  are  the  solutions  of  the  coupled  Jacobi  equations, given  by 


3q{t')  -  m-i  •  Jp(f') 

jp(t')  =  -K{t')3,{t')  (2-7) 


with  the  initial  conditions 


J,(0)  =  0 

Jp(0)  =  I  ,  (2-8) 

where  I  is  the  AT-dimensional  identity  matrix  and  K(t')  is  the  time-dependent  force  constant 
matrix,  K(t')  =  d'^V{t')/dqdq,  evaluated  along  the  stationary  path  V{t')  =  V[qst{,t% 
Clearly,  the  Jacobi  equation  in  Eq.  (2.7)  is  the  same  as  the  equation  of  motion  describing 
an  oscillator  with  a  time-dependent  force  constant  determined  by  the  stationary  trajectory. 

The  nonuniform  semiclassical  formula  in  Eq.  (2.4)  is  valid  as  long  as  the  prefactor  in  Eq. 
(2.6)  remains  finite.  It  happens  at  certain  times  that  two  or  more  paths  may  coalesce  at  a 
focal  point,  or  caustic,  where 
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det  J,(tc)  =  det5q(tc)/5pi  =  0 


(2.9) 


resulting  in  the  divergence  of  the  nonuniform  expression  Eq.  (2.4).  In  that  case,  one  can 
resort  to  more  accurate  uniform  asymptotic  approximations®^  which,  of  course,  assume  a 
more  complicated  form.  Also,  at  a  caustic  tc,  the  number  of  negative  eigenvalues  of  the 
matrix  Jg(fc),  denoted  by  Sign[Jg(tc)],  will  change  depending  on  the  order  of  the  caustic. 
By  keeping  track  of  the  time  evolution  of  Sign[J5(t)],  one  can  express  the  Maslov  index 
explicitly  as 

Kt)  =  y:{SignlJ,K)|  -  SignlJ,K)l}  (2.10) 

k 

where  4  denotes  the  A:-th  caustic  time  as  the  stationary  path  evolves  in  time  from  0  to 
t.  In  fact,  it  can  be  seen  from  Eq.  (2.6)  that  the  Maslov  index  is  simply  the  number  of 
negative  eigenvalues  of  the  second-order  derivative  matrix,  which  will  be  discussed  later  in 
the  context  of  the  initial-value  representation. 

Apart  from  the  difficulties  associated  with  caustics,  the  root  search  problem  poses  a 
formidable  task  in  numerical  implementation  of  Eq.  (2.4).  Unlike  an  initial-value  problem 
where  the  trajectory  follows  a  unique  path  in  phase  space,  the  boundary-value  problem 
requires  one  to  search  for  a  solution  to  Eq.  (2.5)  which  satisfies  both  the  initial  and  final 
conditions  in  Eq.  (2.3),  thus  giving  rise  to  the  possibility  of  multiple  solutions.  For  many- 
body  potentials  there  can  exist  a  very  large  number  of  such  paths  for  longer  time  dynamics. 
One  also  might  obtain  imaginary  paths  in  the  case  of  quantum  tunneling.  The  numeri¬ 
cal  difficulty  associated  with  the  search  for  these  solutions  increases  drastically  with  the 
dimensionality  of  the  system. 

2.  Time  Correlation  Functions 

As  stated  earlier,  many  physical  quantities  of  interest  can  be  related  to  time  correlation 
functions.  In  their  most  general  form,  these  functions  can  be  expressed  as 
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=  j  dqi  j  dq2  j  dqt  j  j  dc{^ 

X  p(q:,q2;/?)(q.|e-^‘/'^|q2)*(q;|e-^‘/'‘|q;)(^^  -  (2-11) 

where  Z  is  the  partition  function,  Z  =  Trexp(— and  p  is  the  canonical  density  matrix 
at  temperature  /?  =  l/ksT,  i.e.,  p(qi,q2;/?)  =  (qile-^^lq2).  While  Eq.  (2.11)  is  a  general 
expression,  for  the  present  discussion  we  will  specialize  it  to  the  case  in  which  the  operators 
A  and  B  are  dependent  on  position  only,  giving 


{i(()  B(0))  = 

Z-> /dq,/dq,/dq.p(qnqn/J)(q.|e-"‘'»lqe)-(q.|e-®''*lq.>/l(q,)B(qi)  ■  (2-12) 

The  case  of  general  operators  depending  on  both  position  and  momenta  will  be  discussed  in 
the  next  subsection  within  the  context  of  the  initial-value  formulation.  By  substituting  the 
semiclassical  formula  in  Eq.  (2.4)  for  the  propagators  into  the  above  expression,  one  obtains 


{A{t)B{0))sc  =  Z-^  I  dqi  I  dq2l  dqtA{qt)B{qi)p{q2,quP) 

1/2 

X  f det|^det|^)  exp(^A5,t/h)  ,  (2.13) 

dpi  dp2j 

where  the  subscript  “sf”  denotes  a  summation  over  both  the  forward  and  backward  station¬ 


ary  paths  and  ASst  =  Sst{qi,qf,t)  -  Sst{q2,qf,t).  In  principle,  time  correlation  functions 
can  be  evaluated  based  on  Eq.  (2.13),  but  such  a  calculation  would  be  fully  vulnerable  to  the 
caustic  and  root  seeirch  problems  described  previously.  Therefore,  the  above  expression  is 
primarily  of  formal  interest.  (It  is  useful,  for  example,  when  one  considers  the  classical  limit, 
cf.  Appendix  A.)  A  much  more  useful  approach  is  based  on  the  initial-value  formulation  of 
semiclassical  dynamics  and  this  will  now  be  discussed. 
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B.  Initial-Value  Formulation 

1.  Propagator 

The  initial-value  representation  of  the  semiclassical  propagator  is  a  recasting  of  the  semi- 
classical  boundary-value  problem  in  terms  of  the  initial  position  and  an  integral  over  the 
initial  momentum.  Since  this  approach  is  formally  equivalent  to  the  Van  Vleck  form,  their 
evaluation  is  formally  equal.  However,  the  initial-value  representation  is  numerically  supe¬ 
rior  since  the  stationary  phase  trajectories  in  the  initial- value  approach  are  determined  from 
initial  momenta  and  coordinates.  The  troublesome  boundary-value  problem  thus  become  an 
initial-value  problem.  Moreover,  the  Van  Vleck  determinant,  which  vanishes  at  the  caustics, 
appears  in  the  numerator,  instead  of  the  denominator,  of  the  semiclassical  expression.  The 
initial-value  representation  is  a  global-time  asymptotic  semiclassical  approximation  which  is 
reducible  to  the  Van  Vleck  formula  by  a  stationary  phase  integration. 

Recently,  Campolieti  and  Brumer^®  presented  an  in-depth  study  of  the  initial-value  for¬ 
malism,  with  an  emphasis  on  a  derivation  of  the  Maslov  indices  and  canonical  transfor¬ 
mations  among  alternative  phase-space  representations.  Their  analysis  follows  a  simple 
procedure  of  concatenating  short-time  propagators  by  sequential  stationary-phase  integra¬ 
tions.  In  the  following  subsection,  we  re-derive  the  initial-value  propagator  from  a  discretized 
perspective  which  leads  to  an  efficient  and  transparent  alternative  for  evaluating  the  Jacobi 
matrices.  These  new  developments  allow  us  to  formulate  an  initial-value  semiclassical  al¬ 
gorithm  for  the  quantum  propagator  in  both  the  adiabatic  and  nonadiabatic  limits.  The 
central  result  of  these  efforts  in  both  cases  is  an  initial-value  expression  for  the  coordinate 
representation  of  the  propagator,  i.e.,®® 

G I  sc  )  ^t'-}  t)  = y  dpi  |det  Jp(t)|  ^  exp['iQ;(qi,  pi,  qt;  ~  i7r/i(t)/2]  ,  (2-14) 

where  the  phase  is  a  canonical  transformation  of  the  classical  action  S,  i.e., 

ot{qi,Pi,qut)  =  Sst[qi,q{t)]t]  +  p{t) -[qt- q{t)]  ,  (2-15) 
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and  the  index  n{t)  is  related  to  the  Maslov  index  by 

^(t)  =  ^(t)  +  Sign[jJ(t)J,(t)]  .  (2.16) 

Here,  the  Jacobi  matrices  3p{t)  and  3q{t)  are  solved  from  Eq.  (2.7)  or  by  the  discrete  approach 
derived  in  the  next  subsection.  The  stationary-phase  condition  determines  the  classical 
trajectory  from  the  usual  initial  conditions  (qi,Pi),  namely, 

q(t)  =  q(qi,Pi;^) 

p(^)  =  p(qi.pi;^) 

which  is  an  initial- value  problem  rather  than  a  boundary- value  problem  as  in  Eq.  (2.3).  Also 
one  sees  from  Eq.  (2.14)  that  a  vanishing  determinant,  det  3p{t),  does  not  lead  to  a  divergent 
prefactor  at  the  caustics. 

Before  proceeding  to  the  next  subsection,  we  note  that  some  care  is  in  order  when  eval¬ 
uating  the  initial-value  propagator  explicitly.  Unlike  the  nonuniform  asymptotic  expression 
in  Eq.  (2.4),  the  initial- value  expression  Eq.  (2.14)  is  nonsymmetric  with  respect  to  the 
exchange  of  the  coordinates  qi  and  q2,  thus  contradicting  the  symmetry  of  the  Green’s 
function  for  a  real  time-independent  Hamiltonian.  To  remedy  this,  one  can  construct  a 
symmetrized  propagator  by  inserting  a  complete  coordinate  basis  set  at  the  half-time,  i.e., 

G(qi,q2;  t)  =  j  dq3(q2le“*^‘^^'^|q3)(q3|e“*^*^^^|qi) 

=  y  dq3G(q2,q3;t/2)G'(qi,q3;t/2)  ,  (2.18) 

where  the  symmetry  property  of  the  Green’s  function  for  time-independent  Hamiltonians 
has  been  used.  In-  the  evaluation  of  time  correlation  functions,  this  symmetrization  is  not 
necessary. 


2.  A  New  Derivation  of  the  Propagator 

The  approach  reviewed  in  the  previous  subsection  involves  solving  the  Jacobi  matrices 
from  the  Jacobi  equation  which,  in  the  case  of  nonadiabatic  dynamics  described  in  Sec.  Ill 
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below,  becomes  a  complex  integro-differential  equation.  To  avoid  this  difficulty,  we  set  out 
to  find  an  alternative  to  evaluate  the  Jacobi  matrices  and  have  thereby  found  it  necessary 
to  derive  the  initial-value  expression  from  a  new  perspective.  Straightforward  and  self- 
contained,  this  derivation  leads  to  a  discretized  expression  for  the  Jacobi  matrices  and  a 
simple  interpretation  of  the  Maslov-like  index.  These  general  expressions  are  applicable  to 
both  adiabatic  and  nonadiabatic  dynamics,  and  are  therefore  in  some  sense  more  general 
than  the  original  expressions. 

To  start,  the  real-time  propagator  is  rewritten  as 


G{qo,qt,t)  =  {qt\e  =  J  dpt{qt\Pt){Pt\e  (2-19) 

where  a  complete  set  of  momentum  state  has  been  inserted.  It  is  then  essential  to  evaluate 
the  position- momentum  propagator  from  the  initial  coordinate  qo  to  the  final  momentum 
Pt,  which  differs  from  the  usual  position-position  propagator  G(qo,  qt;  t)  only  in  the  terminal 
state.  The  position-momentum  propagator  in  the  discretized  path  integral  form  is  given  as 

1  P  /  m  r 

(p.i<=-‘"‘'*iqo> = n  (^)  (2,20) 


where  (j)p  is  a  discretized  canonical  transformation  of  the  action,  given  by 


4>P  = 


2e 


(qi  -  qi_i)  •  m  •  (qi  -  qi_i)  -  e 


-  Pt  ■  qp 


(2.21) 


Here,  P  is  the  discretization  number,  e  is  the  discrete  time  increment  e  =  t/P,  and  m  is,  as 
before,  the  diagonal  particle  mass  matrix. 

The  semiclassical  approximation  is  the  functional  application  of  the  stationary-phase 
approximation.  The  stationary  phase  condition  [dcp/dqilst  =  0  in  the  present  case  determines 
the  discretized  stationary  path  for  i  ^  P  as 

(qi+i  T  qi-i  “  2qi) 


m 


+  VK  =  0 


(2.22) 


and,  for  i  —  P,  &s 

p  =  —  •  (qp  -  qp-i)  -  xVHp  .  (2-23) 

e  2 
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It  is  easy  to  recognize  that  in  the  continuous  limit  Eq.  (2.22)  is  equivalent  to  the  classical 


equation  of  motion,  that  is 


q(i')  +  VE[q(f')]  =  0 


(2.24) 


and  Eq.  (2.23)  imposes  the  terminal  boundary  condition 

Pt  =  P(i)  =  P(qo,Po,t)  ■ 


(2.25) 


Next,  the  quantum  fluctuations  are  evaluated  by  a  second-order  functional  derivative. 


giving 


oqiOqj  e 


(2.26) 


except  for  j  =  P,  which  is  given  by 


aVst  m  e 

-  — - ~^PP  ) 

dc[pdqp  e  2 


(2.27) 


where  K  is  the  time-dependent  force  constant  matrix  Kjj-  —  d^V f dqidqj  evaluated  along 
the  stationary  path. 

The  determinant  of  the  matrix  [d^4>/dq^dqj]  in  the  large  P  limit  can  now  be  defined  as 


_i  8^4) St 

detD„(t)  =  lim  detp  em  — x — 
'  P-oo  dqidqj 


(2.28) 


where  detp  refers  to  the  discretization  of  time  into  P  slices  but  not  to  the  system  dimen 
sionality  N.  This  determinant  is  the  product  of  the  eigenvalues  and  hence  the  phase  /x(t)  of 
Dp(f)  is  determined  by  the  number  of  negative  eigenvalues  of  the  matrix,  i.e.. 


detDp(f)  ^  |detDp(t)|exp[i7r/x(t)]. 


(2.29) 


Thereby,  the  semiclassical  limit  of  the  propagator  in  Eq.  (2.20)  can  be  explicitly  written  as 

(p.|e-‘"‘'''|qo>5c  = - ,  ^  =expliMt)/li  "  (2  30) 

h''^|detD„(()| 

where  (pstit)  =  Sst[qi,q{t)]  t]-p{t)-q{t).  Thus,  the  propagator  in  Eq.  (2.19)  can  be  rewritten 
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G,5c(qo,q.;t)  -  ^/dp„|deU,(()| 

exp{i[(f)st{t) +  p{t)  ■  qt]/h- iTriJ,{t)/2}  ,  (2.31) 

where  a  change  of  variable  from  the  final  to  the  initial  momentum  has  been  carried  out 
which  introduces  the  Jacobi  factor  3p{t).  It  is  proven  in  Appendix  B  that  Dp(f)  is  equal  to 
the  Jacobi  matrix,  i.e., 

Dp(t)  =  3p{t)  .  (2.32) 

This  result  allows  us  to  reach  the  final  expression  of  the  initial- value  representation  given  in 
Eq.  (2.14)  with  pi  in  that  expression  replaced  by  po. 

3.  Time  Correlation  Functions 

By  substituting  Eq.  (2.14)  into  Eq.  (2.11),  it  is  straightforward  to  obtain  an  initial-value 
semiclassical  expression  for  general  time  correlation  functions,  i.e., 

{A{t)B{o))jsc  =  I  I  I  I  I  I  J  dp2  p{qij  q2,  P) 

X  g{Ai,Pi,At]t)9*i<32,P2,qut){<it\A\c^t){q!i\B\qi)  ,  (2.33) 

where  g  is  the  integrand  in  the  initial- value  propagator  in  Eq.  (2.14),  i.e., 

^(qbPbqt;^)  =  |detJp(t)l^'^^exp[lo;(qi,pi,qs;t)/h  —  iTrp.{t)/2].  (2.34) 

In  this  case,  A  and  B  are  general  operators  which  can  depend  on  both  position  and  momen¬ 
tum.  If  one  knows  the  position  matrix  elements  of  these  operators,  and  if  they  are  simple 
products  of  position  and  momentum,  then  one  can  readily  express  the  above  correlation 
function  in  a  more  explicit  form.  For  example,  if  A  and  B  depend  only  on  the  position 
operator,  then  Eq.  (2.33)  simplifies  to  read 

{A{t)B{0))j[sc  —  j  dqij  dq2  j  dqt  j  j  dp2  p{qi,q2,  P) 

X  g{qi,Pi,qt'-,t)9*{q2,P2,qt]t)A{qt)B{qi)  .  (2.35) 
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In  many  cases  of  interest,  however,  the  operators  A  and  B  are  complicated  functions  of 
positions  and  momenta.  In  such  cases,  it  is  better  to  concentrate  one  s  efforts  on  the  time 
correlation  function  written  in  the  form 


{A{t)B{0))  =  ^  j  dqo  j  dqi  j  (iq2  p(qo,  q2;  ^)(q2|^(^)|qi)(qil^(0)|qo)  ,  (2.36) 

where  A{t)  and  5(0)  are  Heisenberg  operators.  The  focus  therefore  shifts  to  deriving  a 
semiclassical  initial-value  expression  for  the  matrix  element  (q2|i(f)|qi).  Through  the  Weyl 
correspondence,  an  operator  can  be  expressed  as®^’®^ 


(q2|i(g,p)lqi)  = -^ydp'^w[(qi +  q2)/2,p']e'^  >  (2.37) 

where  AM.(q,p)  is  the  classical  symbol  corresponding  to  the  quantum  operator  i(q,p).  By 
combining  Eq.  (2.37)  and  the  initial-value  expression  in  Eq.  (2.14),  one  obtains 


(q2|e''^'/M(4p)e'‘'^'/''|qi)/5C  =  j  dq'i  J  dp' Aw  [(q'l  +  q2)/2,  p'] 

X  e^P'-(<-"'^)/"(q;|e-*^‘/"|qi)/5a(q'2|e-‘^‘/^|q2)}5C  ,  (2-38) 


VVe  then  make  a  change  of  variables  to  q  =  (q'l  -  q2)  and  qt  =  (q'l  +  q2)/2,  integrate  over 
the  first  variable  resulting  in  a  delta  function  of  momentum  which  is  also  integrated  over, 
and  finally  arrive  at  the  initial-value  semiclassical  representation  of  the  Heisenberg  operator, 
i.e., 

(q2|e'^‘^^i(q,p)e“'^*'''^|qi)/5C  =  J  dqt  j  dpi  J  dp2Aw  (qt,  [Pi(^)  +  P2(t)]/2) 

X  p(qi,pi,qt;t)5'(q2,P2,qt;t)*  >  (2.39) 


where  the  classical  momentum  in  the  symbol  7lvi/(q,  p)  takes  the  average  value  at  the  end 
of  the  trajectories  which  are  used  in  the  computation  of  p(qi,Piiqt;t)  and  p(q2) P2, qt, i)  • 
Thus,  momentum  symbols  in  the  Weyl  ordering  of  Heisenberg  operators  are  essentially  the 
symmetrized  final  momentum  variables  from  the  initial-value  semiclassical  representation. 
It  can  be  shown  that  Eq.  (2.37)  is  recovered  from  the  t  ^  0+  of  Eq.  (2.39).  With  this 
formulation  of  semiclassical  Heisenberg  operators  in  hand,  the  time  correlation  function  in 
Eq.  (2.36)  can  be  written  in  the  initial-value  representation  as 
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{A{t)B{0))isc=  ^  I  dqo  I  dqy  J  dq2  j  dqt  j  dpi  J  dp2  p{qo,q2]  P) 

X  Aw  iqt,[piit)  +  P2{t)]/‘^)  {qi\B{0)\qo)g{qi, Pi, qt,t)9*{q2,P2,qt,t)  ,  (2.40) 

where  the  zero-time  matrix  element  (qi|5(0)|qo)  has  been  left  in  a  general  form. 

The  form  of  Eq.  (2.40)  also  suggests  a  series  of  approximations  in  which®^ 

/l„(q,p)  =  >l„(q.p)  +  f;i§^/14q,P)  .  P-41) 

n=l 

where  Ao(q,p)  =  Aci{q,p)  is  the  classical  function  of  the  dynamical  variables  q  and  p. 
In  particular,  if  Aw {q,p)  is  a  polynomial  such  a  truncation  will  be  exact  for  sufficiently 
large,  but  finite,  values  of  n.  In  general,  such  a  truncation  is  not  strictly  equivalent  to  the 
semiclassical  approximation  because  the  error  is  in  the  pre-exponential  part.  Nonetheless, 
the  truncation  is  tantamount  to  that  used  by  Wigner®®  in  his  formulation  of  the  leading 
quantum  correction  to  the  classical  partition  function. 

III.  SEMICLASSICAL  THEORY:  NONADIABATIC  DYNAMICS 

Many  advances  have  taken  place  in  the  field  of  nonadiabatic  dynamics  simulation. 

The  theoretical  basis  for  several  algorithms  is  the  Pechukas  theory  of  nonadiabatic  colli¬ 
sions  based  on  the  stationary  phase  approximation  to  Feynman  path  integrals. As  far 
as  the  nuclear  motion  is  concerned,  classical  dynamics  has  been  assumed  in  most  of  the 
nonadiabatic  dynamics  algorithms  based  on  the  Pechukas’  formulation.  Clearly,  neglect  of 
the  quantum  nature  of  the  nuclear  dynamics  is  inadequate  for  treating  light  nuclei  such  as 
protons.  Therefore,  the  combination  of  the  initial-value  semiclassical  approximation  and 
nonadiabatic  dynamics  will  represent  a  more  accurate  description  of  such  systems  and  this 
is  the  focus  of  the  present  section. 

Consider  the  dynamics  of  nuclei  on  a  potential  resulting  from  multiple  electronic  diabatic 
surfaces.  One  then  needs  to  extend  adiabatic  dynamics  to  nonadiabatic  dynamics  to  allow 
for  the  possibility  of  transitions  between  the  different  surfaces.  To  put  the  formalism  in  the 
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most  general  context,  we  consider  the  Hamiltonian  of  a  many-body,  multi-level  system  is 
written  as 

H  =  Ho{q)  +  Hd{q)  ,  (3.1) 

where  q,  as  before,  is  the  collection  of  N  nuclear  degrees  of  freedom  of  the  system  of  interest, 
and  Ho  =  q{t')  ■  m  •  q(t')/2  is  the  kinetic  energy  term  for  the  nuclei.  The  Hamiltonian  Ha 
can  be  explicitly  expressed  in  terms  of  the  elements  ha  (for  the  ith  diabatic  surface)  and  hij 
(for  the  coupling  between  the  fth  and  jth  diabatic  surfaces),  i.e., 

Hd{q)  =  Y^hii  +  Y,J2hij  .  (3.2) 

i  i  j>i 

Here,  the  elements  are  defined  as 

ha  =  |i)  Ki(q)  (^1  ,  (3-3) 

and 

=  \i)V,j{q){j\  +  \j)Vi{q){i\  ,  (3.4) 

where  the  off-diagonal  coupling  elements  satisfy  the  Hermitian  relation  Vij  =  V^*.  The 
potential  energy  terms  Va  in  the  elements  ha  describe  the  diabatic  surfaces,  so  the  above 
formulation  of  the  problem  is  completely  general. 

For  the  Hamiltonian  in  Eq.  (3.1),  the  matrix  element  of  the  nuclear  time  propagator  in 
the  diabatic  basis  reads 

^^.^.(qo.qt;^)  =  J 'Dq{t')exp{iSQ[q{t')]/h}T^u[q{t')]  (3.5) 

where  5'o[q(t0]  the  action  functional  associated  with  the  kinetic  energy  term  Ho  and 
is  the  overlap  between  the  initial  and  final  diabatic  states.  Explicitly,  the  time  evolution 
operator  for  the  diabatic  Hamiltonian  evolves  according  to  the  time-dependent  Schrodinger 
equation 

ih  ^  =  H4q{t')W)  (3.6) 


17 


with  the  initial  condition  u{0)  =  1.  The  transition  amplitudes  are  thus  given  by 


(3.7) 


which  is  a  functional  of  the  nuclear  path  q(t')- 

To  facilitate  the  subsequent  analysis,  the  quantum  average  over  the  diabatic  basis  is 

introduced  as 


(3.8) 


where  the  denominator  is  independent  of  the  variable  t',  and  f{t')  is  in  general  an  operator. 
This  quantum  average  is  carried  out  by  assuming  a  particular  nuclear  path  q{t')  and  is  thus 
a  functional  of  the  nuclear  paths.  Equations  (3.1)-(3.8)  represent  the  exact  formulation  of 
the  nonadiabatic  quantum  dynamics  for  a  given  nuclear  path  q(t'). 

Following  Pechukas’  analysis,  we  apply  the  stationary  phase  approximation  to  Eq.  (3.5) 
and  thus  obtain  the  equation  of  motion  for  the  nuclear  coordinates^^’®^ 

/dH,[q{t')]\ 


m  ■  q(f')  =  - 


(3.9) 


\  5q(i0  / d  ’ 

which  is  to  be  solved  together  with  Eqs.  (3.6)  — (3.8)  to  obtain  the  nonadiabatic  stationary 
solution(s).  Clearly,  the  coupling  between  the  diabatic  state  propagation  and  the  station¬ 
ary  trajectory  imposes  the  self-consistency  condition  for  the  solution  of  the  nonadiabatic 
trajectories. 

With  the  stationary  solutions  in  hand,  the  semiclassical  expression  for  the  propagation 
is  given  as 


(3.10) 


which  is  the  same  as  the  adiabatic  semiclassical  expression  in  Eq.  (2.14)  except  that  the  Sst 
is  now  defined  as 


Sf,i.,st{qi,fit]t)  =  5'o,st(qi,qt;i)  -  ih\nT^„[qst{t)]- 


(3.11) 
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so  that  Mt)  =  S,.,st  -  p{t)  •  q{t)  and  a{t)  =  +  p(t)  •  [q*  -  q{t)].  The  prefactor 

Jp(t)  obeys  the  Jacobi  equation  which  can  be  obtained  from  taking  the  partial  derivative  of 
Eq.  (3.9)  with  respect  to  the  initial  values  po  or  qo,  giving  the  Jacobi  matrix  equation 


m- J,(0  +  (V'- V'i^d)dJ,(0 

_  1  f' dt”[{{VHd)u{t',t"){V'Hd))d  -  {VHd)d  ■  {V'Hd)d]Jg{t”)  =  0  ,  (3.12) 

h  Jo 

where  V'  =  d/dq{t'),  V  =  dfdq{t"),  and  m  •  j,  =  Jp.  The  latter  equations  are  the 
nonadiabatic  analog  to  the  adiabatic  Jacobi  equations  of  Eqs.  (B6)— (B7)  [cf.  Appendix  B]. 

Solving  the  integral-differential  equation  in  Eq.  (3.12)  is  a  formidable  task.  In  addition, 
the  fact  that  the  Jacobi  matrices  are  generally  complex  introduces  ambiguities  in  defining 
the  Maslov  index  n{t).  To  circumvent  these  difficulties,  we  resort  to  a  discretized  expression 
for  the  Jacobi  matrix  which  is  equivalent  to  solving  the  Jacobi  equation.  To  be  more  explicit, 
the  Jacobi  determinant,  detJp(t),  can  be  evaluated  in  a  discretized  format  as  in  the  previous 


section,  giving  in  the  nonadiabatic  case 


St 


=  ^  -  Sij+i  -  5,j_i)  -  e 


dqidqj 


'  d^Hd[q{t')]\ 

.  5q?  L 


(3.13) 


except  for  f  =  j  =  P,  which  is  given  by 

ay,,  _  m  .  (3.14) 

dqpdqp  t  2  \  5qp  2 

The  quantum  fluctuation  correlation  matrix  Cij  here  is  given  by 

/SH4q(t')]  ,  ,  dHi[<i{t')]\  /dHilq(t')]\  /  3gj[q(t')]  \ 

\  Sq.  /A  aq,  // 

The  dimensionality  implicit  in  the  above  equations  is  such  that  d'^(f>st/dqidqj  is  a  matrix  of 
dimension  N  x  P.  After  taking  the  limit  P  ^  oo,  we  obtain  the  explicit  expression  for  the 
prefactor 


detJp(t)  =  lim  det  em' 

P— )-oo 


dqidqj 


(3.16) 


Here,  the  determinant  denotes  the  product  of  eigenvalues  which  are  complex  and  defined 
counter-clockwise  in  the  complex  plane.  Therefore,  the  Maslov-like  index  /.t(t)  equals  the 
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sumniation  of  the  phase  angles  of  the  complex  eigenvalues  of  the  discretized  second  order 
derivative  matrix  in  Eqs.  (3.13)  and  (3,14)  in  the  large  P  limit.  The  above  derivation 
represents  a  different  approach  for  calculating  the  Jacobi  matrices  and  the  Maslov-like  index 
without  solving  the  Jacobi  equations.  It  is  applicable  in  both  the  adiabatic  and  nonadiabatic 
limits. 


IV.  NUMERICAL  ALGORITHMS 

Though  the  initial- value  semiclassical  expression  in  Eq.  (2.14)  represents  a  significant 
simplification  of  the  exact  path  integral,  the  integrand  g  of  Eq.  (2.34)  as  a  function  of  the 
initial  momenta  is  oscillatory  and  thus  does  not  render  itself  to  simple  integration  schemes. 
Since  the  usual  Monte  Carlo  method  is  not  applicable  to  integrate  such  complex  expo¬ 
nents,  one  has  to  introduce  a  positive  definite  weight  function  such  that  the  integration 
domains  which  dominate  the  integral  will  be  sampled  preferentially  over  those  which  barely 
contribute  because  of  phase  cancellations.  Indeed,  such  an  approach,  termed  stationary 
phase  Monte  Carlo,  has  been  proposed  and  applied  in  some  quantum  dynamical  path  inte¬ 
gral  simulations. Here,  we  will  describe  a  simplified  version  of  the  stationary  phase 
method  as  it  applies  to  the  present  semiclassical  formalism. 

In  general,  consider  a  one-dimensional  integral  of  the  form 

1(h)  =  r  dx  ,  (4.1) 

J — oo 

which  is  a  generic  integral  having  a  complex  exponent.  If  fi  is  small,  the  integral  is  dominated 
by  regions  where  the  phase  (j){x)  is  stationary,  i.e.,  where  =  0.  In  the  non-stationary 
regions,  the  complex  exponential  is  highly  oscillatory,  thus  effectively  canceling  the  contribu¬ 
tion  from  those  regions.  Therefore,  it  is  advantageous  to  introduce  a  weight  function  which 
suppresses  the  oscillatory  integrand  and  thus  acts  effectively  as  a  filter.  A  simple  example 
is  a  Gaussian  filter,  defined  as 

Wfix)  =  exp{-t[fi\x)]^}  ,  (4-2) 
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where  e  is  the  filter  parameter.  Then,  the  integration  in  Eq.  (4.1)  becomes 

/CC 

dx  ,  (4.3) 

'CO 

and,  for  example,  an  expectation  value  is  approximated  as 

A  general  form  of  filters  and  the  relationship  between  the  filtered  integral  in  Eq.  (4.3)  and  the 
exact  one  in  Eq.  (4.1)  have  been  analyzed  in  detail  by  others. It  should  be  noted  that 
there  is  some  flexibility  in  the  definition  of  the  pre-exponential  factor  of  the  filter  function 
depending  on  the  final  numerical  target,  but  these  factors  will  cancel  in  the  calculation  of 
expectation  values  or  time  correlation  functions  as  outlined  below. 

Although  approximations  to  the  exact  integrals,  Eqs.  (4.3)  and  (4.4)  converge  much 
faster  than  the  original  integrals.  If  the  parameter  t  is  small,  the  filtered  integral  is  closer 
to  the  exact  one  but  takes  much  longer  to  converge;  in  the  case  of  e  0,  the  exact  integral 
is  recovered.  If  the  parameter  e  is  large,  the  filtered  integral  is  localized  near  stationary 
points  and  thus  ignores  fluctuations  away  from  those  points.  A  proper  choice  of  the  filter 
parameter  is  indeed  crucial  for  carrying  out  an  accurate  and  efficient  calculation.  Given  a 
required  level  of  convergence,  the  filtered  integral  exhibits  poor  statistics  for  e  <  Cmin:  but 
it  may  deviates  substantially  from  the  exact  value  for  e  >  Cmax-  Thus,  the  optimal  choice  of 
e  is  located  in  the  intermediate  region,  Cmin  <  e  <  ^max,  where  the  filtered  integral  becomes 
both  stable  and  accurate. 

Following  the  above  discussion  of  the  stationary-phase  filter  method,  it  is  now  specialized 
to  treat  the  integrations  in  the  initial-value  representation  of  semiclassical  time  correlation 
functions  explicit  in  Eq.  (2.33).  Since  the  general  term  Q!(qi, Pi, qp,  t)  defined  by  Eq.  (2.15) 
is  the  phase,  the  filtered  propagator  is  given  by 

{q2\e~"^''^^\qi)isc  =  ^  j  dpig{qi,Pi,q2-,t)W,{pi)  ,  (4.5) 

where  the  filter  can  be  defined  as 
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W.p,  =  exp 


{da/dpiY  =  exp  ep{Jp{t)[qt  —  Q{i)]V) 


(4.6) 


and  a  one-dimensional  notation  has  been  adopted  here  for  notational  simplicity.  It  is  also 
straightforward  to  write  down  the  filtered  expression  for  the  semiclassical  time  correlation 
function.  For  example,  the  case  of  position-dependent  operators  is  given  from  Eq.  (2.35)  by 


{A{t)B{0))isc  = 

/  dqi  f  dq2  f  dqt  J  dpi  f  dp2  p{qi,q2',  P)  Kg{t)'^e,pi^e,P2^i,qf^{Qt) B {qi)  ^ 

/  dqi  J  dq2  f  dqt  f  dpi  f  dp2  p{qi,q2]  P)  Kg{t)^Ve.pi^€,p2We,qt 

where  Kg{t)  g{qi, Pi,  qt]  t)*g{q2,P2,  qt,t),  and  a  filter  is  also  applied  here  to  the  qt  integra¬ 
tion,  i.e, 

=exp{-eg[5(ai -02)/%]^}  =exp{-eg[pi(f) -P2(i)]^}  •  (4-8) 

In  Eq.  (4.7) ,  the  quantum  density  matrix  element  p{qi,q2',  P)  also  provides  a  natural  filter 
for  the  integration  over  the  variables  qi  and  q2-  The  functions  A{qt)  and  B{q\)  in  Eq.  (4.7) 
may  also  aid  in  the  filtering,  depending  on  their  form.  Note  that  the  denominator  of  Eq. 
(4.7)  equals  the  partition  function  Z,  but  it  has  been  written  so  that  the  overall  equation  is 
amenable  to  a  Monte  Carlo  algorithm.  It  should  also  noted  that  the  trajectories  q{t)  in  Eq. 
(4.6)  which  contribute  to  and  VK,p,  in  Eq.  (4.7)  are  different  (i.e.,  they  have  different 
initial  momenta)  and  therefore  must  be  treated  as  such. 


V.  NUMERICAL  EXAMPLES 
A.  Propagator  for  a  Solvable  Potential 

To  demonstrate  the  feasibility  and  accuracy  of  the  initial-value  semiclassical  approxi¬ 
mation  and  the  stationary-phase  filter  technique,  we  first  present  a  numerical  study  of  a 
solvable  potential,  given  by 

m  =  ^ 
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with  m  =  1.0  and  h  =  1.0.  The  real-time  propagator  of  this  potential  is  given  in  a  closed 
form  by® 


(5.2) 


where  Ii,  is  the  modified  Bessel  function  of  index  u  —  3/2.  In  Fig.  1  is  plotted  the  unfiltered 
integrand  g  of  Eq.  (2.34),  the  filtered  integrand  in  Eq.  (4.5),  and  the  filter  in  Eq.  (4.6)  as  a 
function  of  the  initial  momentum  for  q,  =  3,qt  =  i,t  =  3.  The  filter  parameter  in  Eq.  (4.6) 
was  taken  as  =  0.01.  It  can  be  seen  clearly  from  the  figure  that  the  integrand  is  highly 
oscillatory  except  near  the  origin.  The  filter  selects  two  stationary  regions,  one  to  the  left 
of  the  origin,  the  other  to  the  right  of  the  origin,  which  indeed  correspond  to  two  possible 
classical  trajectories,  the  direct  path  and  the  indirect  path  bounced  from  the  repulsive  wall. 

The  squared  amplitude  of  the  time  propagator  for  the  potential  in  Eq.  (5.1)  with  the 
same  parameters  as  in  Fig.  1  is  plotted  as  a  function  of  time  in  Fig.  2,  where  the  exact  and 
semiclassical  results  are  represented  by  a  solid  and  dashed  lines,  respectively.  The  initial 
momentum  was  integrated  on  a  200-point  grid  from  —10  to  10  with  a  filter  parameter  of  Cp  — 
0.01.  Despite  the  small  discrepancies  due  to  the  nature  of  the  semiclassical  approximation, 
good  agreement  with  the  exact  result  is  achieved. 


B.  Anharmonic  Quantum  Oscillator:  Position  Correlation  Function 

In  this  subsection,  the  initial-value  semiclassical  method  is  used  to  compute  the  position 
correlation  function  for  an  anharmonic  potential,  defined  as 

with  m  =  1.0,  h  =  1.0  and  /3  =  1.0.  The  position  correlation  function  was  computed  with 
filters  on  the  momentum  integrations  having  a  value  of  Cp  =  0.1.  The  thermal  density  matrix 
was  calculated  by  the  numerical  multiplication  method  (NMM),^^  and  the  coordinates  and 
momenta  were  integrated  on  grids.  The  numerically  exact  correlation  function  was  obtained 
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from  a  harmonic  oscillator  basis  set  calculation.  In  Fig.  3,  the  real  parts  of  the  exact  and 
semiclassical  time  correlation  functions  are  plotted  as  functions  of  time.  The  semiclassical 
approximation  is  more  accurate  for  time  correlation  functions  than  for  propagators,  probably 
because  correlation  fnnctions  result  from  a  thermal  average  of  forward  and  backward  real¬ 
time  propagators  and  are  thereby  less  sensitive  to  errors  introduced  by  the  semiclassical 
approximation.  In  general,  the  agreement  between  the  exact  and  semiclassical  results  is 
excellent  for  this  system. 


C.  Double  Well:  Reactive  Flux 

As  a  final  example,  the  flux-flux  correlation  function  was  calculated  for  a  double  well 
potential,  defined  as 

n<l)  =  -lf  +  (5-4) 

with  m  =  1.0,  p  =  1.0  and  h  =  1.0.  The  quantum  dynamics  of  a  double  well  exhibits 
coherence  at  low  temperature,  thermal  activation  at  high  temperature,  and  assumes  an 
exponential  decay  in  the  presence  of  dissipation.  Miller  and  co-workers  have  shown  that 
thermal  rate  constants  can  be  obtained  from  the  time  integration  of  the  flux-flux  correlation 
function,  defined  as^^ 

CpFit)  =  Tr  ,  (5-5) 

where  the  flux  operator  is  given  by 

F  = -^\p6{q  -  Qb)  +  6{q  -  qb)p]  (5-6) 

2m 

with  qb  defined  as  the  position  of  the  dividing  surface.  To  be  more  explicit,  we  define  a 
complex  time  propagator  as 

=  (5-7) 

such  that  the  flux-flux  correlation  function  can  be  expressed  as^^ 
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Cpp{t)  =  -  lim  \\pG{qt,q]f3,t)\^  +  'ReppbG{qb,q]f3,t)G*{qb,q-,P,t)]  (5.8) 

Q-^Qb  2m^  ^ 


where  pb  is  the  momentum  operator  acting  on  qb-  Although  the  flux-flux  correlation  function 
depends  upon  the  choice  of  the  dividing  surface,  the  rate  constant  is  a  physical  quantity 
independent  of  qb-  For  convenience,  we  choose  =  0  so  that  the  first  term  in  Eq.  (5.8) 
vanishes  for  a  symmetric  barrier.  In  Fig.  4,  the  semiclassical  value  of  Cprit)  is  plotted  for 
the  double  well  and  compared  with  the  exact  result  obtained  from  a  basis  set  calculation. 
Again,  excellent  agreement  is  obtained. 


VI.  CONCLUDING  REMARKS 

In  this  paper,  we  have  discussed  the  semiclassical  formulation  of  quantum  dynamical 
time  correlation  functions  and  have  also  investigated  the  numerical  feasibility  of  the  semi¬ 
classical  approximation  for  calculating  such  functions.  We  demonstrated  the  reduction  from 
the  exact  quantum  time  correlation  function  to  a  nonuniform  boundary-value  semiclassical 
expression,  then  to  a  global-time  initial-value  semiclassical  representation,  and  finally  to  the 
limit  of  electronically  nonadiabatic  quantum  nuclear  dynamics.  Much  of  this  formulation 
was  accomplished  with  the  help  of  a  discrete  approach.  The  resulting  discrete  initial-value 
representation  of  the  semiclassical  approximation  proves  to  be  advantageous  for  the  im¬ 
plementation  of  semiclassical  dynamics  since  the  global-time  formula  avoids  the  problems 
associated  with  caustics  and  root  searches. 

The  studies  presented  in  this  paper  are  not  only  instructive  and  revealing,  but  they 
also  serve  as  the  formal  basis  for  numerical  algorithms.  To  achieve  numerical  efficiency 
in  such  algorithms,  a  stationary-phase  filter  technique  was  introduced  into  the  method  to 
effectively  suppress  the  oscillatory  region  of  the  integrations.  Several  examples  were  tested 
with  the  semiclassical  method  and  compared  with  the  exact  results  obtained  from  basis- 
set  calculations.  These  studies  clearly  demonstrated  the  feasibility  and  accuracy  of  the 
algorithm. 
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With  the  results  of  the  present  work  in  hand,  non-trivial  applications  of  the  semiclassical 
theory  should  be  within  reach.  For  example,  the  combined  use  of  the  initial-value  semiclas¬ 
sical  formalism  and  the  stationary  phase  Monte  Carlo  technique  should  allow  us  to  calculate 
time  correlation  functions  for  realistic  many-body  systems,  particularly  for  one  (or  a  few) 
quantum  particles  in  a  classical-like  environment.  It  will  also  be  very  interesting  to  apply 
this  algorithm  when  the  semiclassical  nuclear  dynamics  of  such  systems  must  be  treated 
nonadiabatically.  For  more  complex  quantum  systems  in  the  condensed  phase,  the  semiclas¬ 
sical  method  can  also  be  used  to  accurately  propagate  an  “important”  quantum  subsystem 
(e.g.,  solute)  under  the  influence  of  an  approximate  quantum  environmental  force  calculated 
by  or  from  a  quantum  Gaussian  bath.  We  have,  in  fact,  developed  a  simple  and 

flexible  scheme^^  to  generate  the  quantum  forces  in  the  latter  scenario  which  can  readily  be 
incorporated  into  the  semiclassical  methodology  described  in  this  paper.  These  and  other 
developments  should  greatly  facilitate  our  ongoing  efforts  to  numerically  simulate  a  wide 
range  of  complex  quantum  dynamical  processes  in  condensed  matter. 
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APPENDIX  A:  THE  CLASSICAL  LIMIT  OF  QUANTUM  TIME  CORRELATION 

FUNCTIONS 

To  reveal  the  relationship  between  classical  and  quantum  time  correlation  functions,  one 
can  introduce  into  Eq.  (2.13)  a  set  of  collective  coordinates^'* 

q(t')  =  qf{t')  -  qb{t') 

q{t')  =  lq,{t')  +  qmft 

with  corresponding  collective  momenta 

p(t')  =  -  pb{t') 

p(^0  [p/(^0  +  P6(i0]/2  ) 

where  q/it')  and  qb{t')  are  the  forward  and  backward  classical  paths,  respectively,  and  Pf{t') 
and  Pb{t')  are  the  corresponding  momenta.  Note  that  in  the  context  of  the  time  correlation 
function  given  by  Eq.  (2.13)  the  initial  collective  coordinates  are  given  by  q(0)  =  (qi  +  q2)/2 
and  q(0)  =  qi  —  q2  aiid  the  final  collective  coordinates  are  given  by  q{t)  =  qt  and  q{t)  =  0. 
Assuming  that  the  path  difference  q{t')  is  small,  one  can  expand  the  action  difference  ASst 
in  Eq.  (2.13)  to  linear  order  in  q  such  that 

ASst  =  Sst{quqt;t)  -  Sst{q2,qut)  ^  -Po  •  qo  (^3) 

where  po  =  Pst(O)  and  qo  =  q5t(0).  Furthermore,  the  difference  of  the  forward  and  backward 
stationary  paths  can  be  ignored  in  the  non-exponential  factor  of  Eq.  (2.13)  and  a  Jacobian 
transformation  can  be  performed,  giving 

j  dqt  |det(9qt/(9por^  =  j  dpo  (-^.4) 

which  changes  the  semiclassical  boundary- value  problem  to  an  initial-value  problem.  Putting 
all  the  pieces  together  and  omitting  the  irrelevant  indices,  we  have 

{A{t)BiO))w  =  /^q/  dpW{q,p-J)A[qait)]B[qam  ,  (A5) 
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where  qci{t)  is  the  classical  trajectory  and  W{q,p)  is  the  well-known  Wigner  distribution 
function,  defined  as®^ 

VF(q,  p]P)  =  ^  I  dq  p(q  +  q/2,  q  -  q/2;  p)  .  (A6) 

Assuming  the  high^temperuture  approximation  of  the  canonical  density  matrix,  i.e., 

Mq..®;ffl=(^)  'exp{-^(q.-q,)-m.(q,-q,)-/JV^(5i±^)).  (A7) 

the  Wigner  distribution  function  (A6)  reduces  to  the  Boltzmann  distribution  function.  The 
classical  time  correlation  function  is  then  completely  recovered,  giving 

{A{t)B{0))ci  =  ^  [-PHci{p,<i)]  A[qci{t)]B[qci{0)]  (A8) 

with  H{p,  q)  and  Zd  being  the  classical  Hamiltonian  and  partition  function,  respectively. 
Thus,  classical  dynamics  results  from  the  high  temperature  approximation  of  the  density 
matrix  and  the  linear  expansion  of  the  action  difference  of  the  forward  and  backward  sta¬ 
tionary  paths  in  the  semiclassical  expression  for  the  time  correlation  function.  The  reduction 
to  the  classical  limit  starting  from  the  initial-value  approximation  does  not  differ  from  the 
above  derivation.  The  Weyl  ordering  of  operators  simplifies  the  analysis,  though  the  classical 
limit  does  not  depend  on  operator  orderings. 

Clearly,  Eq.  (A5)  implies  that  a  quasiclassical  dynamics  can  be  constructed  based  on 
the  Wigner  distribution  function.  On  the  other  hand,  the  fact  that  the  Wigner  distribution 
function  is  not  positive-definite  complicates  the  interpretation  of  the  Wigner  distribution  as 
a  phase-space  quantum  distribution  function.  However,  one  might  adopt  a  coarse-grained 
Wigner  distribution,  such  as  the  Husimi  distribution  function,'"’  as  a  quantum  analogy  to 
the  classical  Boltzmann  distribution  function. 
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APPENDIX  B:  DISCRETE  DERIVATION  OF  THE  JACOBI  EQUATION 

For  convenience,  the  derivation  in  this  appendix  is  presented  for  one  degree-of-freedom 
as  the  multidimensional  generalization  presents  no  special  difficulties.  To  start,  we  define 
the  determinant  of  the  following  two  matrices  as 


2-e^Wi  -1  0 

-1  2-£^W2  -1 


Dq{P)  =  det 


-1  2-^Wp-2  -1 

0  -1  2  - 


which  appears  in  the  discrete  expression  for  the  position-position  propagator,  and 


Dp{P)  =  det 


2-e^Wi  -1  0 

-1  2-e^W2  -1 


■1  2-e^Wp-i 


■1  1  - 


which  appears  in  the  position-momentum  propagator.  Note  that  the  determinant  here  is 
the  product  of  the  eigenvalues  and  thereby  does  not  imply  the  absolute  value.  It  is  then 
easy  to  observe  the  following  iterative  relations 

Dg(P  +  1)  =  (2  -  €^Wp)Dg{P)  -  Dg{P  -  1)  (B3) 


Dp{P)  =  fl  -  {ey2)Wp\  Dg{P)  -  Dg{P  -  1) 


Combining  Eqs.  (B3)  and  (B4),  we  have 


Dp{P)  =  -[Dg{P+l)-Dg{P-l)] 
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To  associate  the  above  difference  equation  to  differential  equations,  we  take  the  limit  e  0, 
define  the  continuous  variable  t  =  eF,  and  introduce  continuous  variables  J,{t)  =  eDg{P)/m 
and  Jp{t)  =  Dp{P).  With  these  definitions  in  hand,  we  can  rewrite  Eqs.  (B3)  and  (B5)  as 

Jg  +  W{t)Jg  =  0  (B6) 

and 

mjq  =  Jp 

with  the  boundary  conditions  specified  as  J,(0)  =  0  and  Jp{0)  =  1.  Obviously,  Eqs.  (B6)  and 
(B7)  are  identical  to  the  Jacobi  equations  if  W{t)  is  specified  as  mW{t)  =  (fV[qst{t)]/dq  , 
therefore,  J,  and  Jp  are  the  usual  Jacobi  variables,  i.e.,  Jqit)  =  dq{t)/dpi  and  Jp{t)  = 
dp{t)/dpi. 
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FIGURES 

FIG.  1.  The  unfiltered  integrand  g  in  the  initial-value  semiclassical  approximation  (dashed 
line),  the  filtered  integrand  in  Eq.  (4.5)  (solid  line),  and  the  filter  function  (bold  line)  plotted  as 
functions  of  the  initial  momentum  for  the  l/q^  potential  with  qi  =  3,  92  =  4,  and  t  =  3. 

FIG.  2.  A  plot  of  |G/5c(3,4;t)p  versus  time  (dashed  line)  for  the  l/q^  potential.  The  exact 
result  is  shown  for  comparison  by  the  solid  line. 

FIG.  3.  The  real  part  of  the  initial-value  semiclassical  position  correlation  function  (solid  di¬ 
amonds)  for  the  potential  in  Eq.  (5.3)  at  a  temperature  /?  =  1.0.  The  numerically  exact  result 
obtained  from  a  basis  set  calculation  is  shown  by  the  solid  line. 

FIG.  4.  The  initial-value  semiclassical  flux-flux  correlation  function  (solid  diamonds)  for  the 
potential  in  Eq.  (5.4)  plotted  along  with  the  numerically  exact  result  (solid  line)  obtained  from  a 
basis  set  calculation 
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